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Abstract. In this paper we propose an algorithm for the numerical 
solution of arbitrary differential equations of fractional order. The algo- 
rithm is obtained by using the following decomposition of the differential 
equation into a system of differential equation of integer order connected 
with inverse forms of Abel-integral equations. The algorithm is used for 
solution of the linear and non-linear equations. 



1 Introduction 

In opposite to differential equations of integer order, in which derivatives depend 
only on the local behaviour of the function, fractional differential equations ac- 
cumulate the whole information of the function in a weighted form. This is so 
called memory effect and has many applications in physics p0[ , chemistry ||, 
engineering ||, etc. For that reason we need a method for solving such equa- 
tions which will be effective, easy-to-use and applied for the equations in general 
form. However, known methods used for solution of the equations have more 
disadvantages. Analytical methods, described in detail in J7|jl2|, do not work in 
the case of arbitrary real order. Another analytical method which uses the 
multivariate Mittag-Leffler function and generalizes the previous results, can be 
used only for linear type of equations. On the other hand, for specific differential 
equations with oscillating and periodic solution there are some specific numerical 
methods H[D|[L4j . Other numerical methods allow solution of the equations 
of arbitrary real order but they work properly only for relatively simple form of 
fractional equations. 

Let us consider an initial value problem for the fractional differential equation 

ai T D^y(t) + a 2 T D? 2 y{t) +... + a n T D?"y(t) + a n+1 y(t) = f(t) (1) 

connected with initial conditions 

y (k) (T) = b k , (2) 

where < t < T < oo, a, S H, a\ ^ 0, 014 £ H+, i = 1, . . . , n, rrii — 1 < < nij, 
mi G IN, ai > ai+i for I = 1 , . . . , n — 1, b k G H, k = , . . . , mi — 1, f(t) is a given 



function defined on the interval [0, T], y(t) is the unknown function which is the 
solution of eqn. (|l|). 

The fractional derivative operator D ai is defined in the Ricmann-Liouvillc 
sense [fL2| . We also have other definitions of the operator like Caputo ||, 
Grunwald-Letnikov ||, and Weyl-Marchaud |fl2|| . Regarding to the Riemann- 
Liouville operator we can define it as the left-side operator 

( D « v)(t) l £L f ym (3) 

U * V)[ ' ~ r(m - a) df» J (t- ' [6) 

a 

where a < t < T < oo, m— \ < a < m, m € IN. In our study, we also use integral 
operators defined in the Riemann-Liouville sense Jl^]. When y(t) 6 Li(a, b) and 
a > then left-side integral operator is defined as 

t 

a 

More informations concerning the operators properties one can find in litera- 
ture [j7p|]l2|| . Applying definitions @ and (||) we have the following property 

( a D?y)(t)=D m { a ir- a y)(t), (5) 

where D m represents typical derivative operator of integer order m. According 
to H in the fractional integrals the composition rule occurs 

„Wf !/(«)) = al? +f) y{t) = aI? +a V (i) = al?{al?y{t)) ■ (6) 

In the general case the Riemann-Liouville fractional derivatives do not commute 

( a D?( a D?y))(t) * ( a D? +0 )(t) + ( a D?( a D?y))(t) . (7) 

Extending our considerations, the integer operator commutes with the fractional 
operator 

(D?( a D?y)){t) = ( a D? + °<)(t), (8) 

but the opposite property is impossible. The mixed operators (derivatives and 
integrals) commute only in the following way 

{ a Bn a lly)){t) = { a D?-P){t), (9) 

but they do not commute in the opposite way. 

We turn our attention to the composition rule in two following ways. First of 
all we found in literature Jl^,|ll| the fact, that authors neglect the general prop- 
erty of fractional derivatives given by eqn. (^). Regarding to solution of fractional 
differential equations we will apply above properties in the next sections. 



2 Mathematical background 



In this section we concentrate on describing the method which is decomposition 
of fractional differential equation into a system of one ordinary differential equa- 
tion and inverse forms of Abel-integral equations. We focus on decomposition 
ways of arbitrary fractional differential equation ([!]) with initial conditions (|^) . In 
our considerations we classify the fractional differential equations in dependence 
on the fractional derivative occurrence as: 

— one-term equations, in which D a occurs only one, 

— multi-term equations, in which D a occurs much more times. 



Corollary 1. An arbitrary one-term equation 

a x T D?y{t) + a 2 y{t) = f(t) , 0<m-l<a<m, melN, (10) 

in which ax ^ and initial conditions y^ k '(r) — bk, (k = , . . . , m — 1) should 
be taken into account, can decomposed into the following system 



'D?z(t) = -%y(t) + ±f(t) 

y {0) {t)^Y:^b k ^ + rDT' a z{t) 
y (1) (t) = + T Dr- a+1 z(t) 

It-rY 



(11) 



/"'-!) (<) =b 



ra — X 



(m-1)! 



r-LJ + 



Proof. We use the property ([5]) and we introduce a new variable z(t) which we 
call a temporary function. For such introduction we have 

T D?y(t) = Dl n z(t) and z(t) = (r^v)®- 

The new variable represents the Abel integral equation which solution, found in 
literature is the left inverse operator 

y(t) = ( T D?z)(t). 

In dependence on the integer order m we introduce initial conditions bk which 
are multiplied by a term k , ■ The term issues from a kernel of the left-side 
Riemman-Liouville operator (||). □ 



In the system of eqns. (|1 1|) , the first equation is the differential equation of m 
integer order and the next equations represent inverse forms of Abel-integral 
equations. We also found in literature HQ different or similar approaches for 
solution of the eqn. (|To|). 

In mathematical point of view, more interesting is the second class called 
multi-term equations. Let us consider the eqn. (|l|) with initial conditions (||) . As 
in previous class we use the rule (|5|) . Following the previous corollary we apply 



■ ai Dr 


21 (i) + 






+a n - 


ny(t) = f(t) 


Zl(t) = 


(rir 


~ ai y)(t) 


z 2 (t) = 


(rir 


- a2 y)(t) 


, z n(t) = 







a system of new variables Zi{t) — ( T Il ni ai y){t) and we decompose eqn. (|l|) into 
a system of the following equations 



(12) 



Such system is not a final form and requires additional transformations. In multi- 
term class we distinguish two sub-classes depending on oti or follows the decom- 
position on m,: 

— independent subclass, in which rm are different derivatives orders of integer 
type mi ^ m 2 :+i for i = 1, . . . ,n — 1, 

— dependent subclass, in which one can find some dependence between the 
derivatives order of integer type, e.g. m\ = m 2 = rn 3 or m\ — m±, etc. 

Proposition 1. Taking into consideration the independent subclass it is possible 
to formulate the unique solution of eqn. ( fli[ ) in the following form 

■ D^z x {t) = -± (03 D™z 2 {t) + ... + a n £>?*"*»(*)) + 
+ ± (-a n+1 y(t) + f(t)) 
z 2 (t) = D r t ni - ,n2 - ai+a2 Zl (t) 

Z3 (t) = Dr- m3 - ai+aa Zl (t) 



Zn(t) = oDT 1 



-m n —a.± -\-a r . 
(t~r) k 



y (1) (t) = r^ 1 h^. 



*l(t) 

+ oA mi " ai ^i(*) 



(13) 



mi-ai + l 



Z X {t) 



k V \ l ) — 0»ni-l (mi-l)! + a±J t 



2mi-ai-l 



Zi(t) 



Proof. In this proof we try to show a dependence of temporary functions Zi (i) 
(i = 2, n) on the temporary function z\(t). We assume, that initial conditions 
bk = 0. In such way we can easily construct the proof. We observe that the Abel 
integral operators in ( |l2| ) have the order < mi — 014 < 1. We consider the 
inverse form of such operators and for two operators we have 

T D^z 2 (t)= T D] Zl (t), (14) 

where < /3 < 1 and < 7 < 1 respectively. Using the property (JsJ) we obtain 

Neglecting the operator D\ we can obtain the following formula 



rll^Z^t) 



1-7 



zi(t) + Ct 



(15) 



where G\ 2 is an arbitrary constant. Additionally, we multiply the operator T lf 
for eqn. ( |T5| ) and applying the property (j(^) we obtain 

r l\z 2 {t) = T I^ +f3 Zl (t) + T jf Ci l2 . (16) 

Differentiating both sides by the operator D\ and applying the property (|^) (see 
literature we obtain the final dependence 

z 2 (t) = T D7^ Zl {t) + Cx, a • r(1 + ^. fl+g - (17) 

In general way we state that Zi(t) [i = 2, ...,n) depends on z±(t) together with 
additional term represented by the constants Ci j. □ 



Proposition 2. All constants Cu = (i = 2, ...,n) m egn. (|i_/[) are equaled to 
zero for r = 0. 

Proof. Taking into consideration the formula ( |l7| ) we can calculate the constants 

r(i + /?)-t 1 +/ 3 / fl n 

Ci,* = - [Zi{t) - T D{ ' z^t) J . 

Additionally, putting the initial condition t = t = we obtain all C\j — 0. □ 

Remark 1. According to previous considerations, the initial conditions for tem- 
porary functions Zi(t) satisfy the following dependence 

«i(r = 0) = 0, i = l,...,n. (18) 

In practical point of view given by formula (|l^), we can solve ordinary differential 
equation of integer type with variable coefficients under zeros initial conditions 
of temporary functions z\ (t) . 

Let us consider the next subclass of a fractional differential equation. We 
assume, that some integer orders rm are the same. This may happen in the 
case when some ai and ai+i, (i = 1, ... ,n — 1) belong to the same integer 
values m, = JTij+i. We introduce a parameter r which denotes a number of 
temporary functions having the same derivative operator D mi of integer type. 
If mi — m 2 = . . . = m r for different temporary functions Zj(t), [i = 1, . . . , r) 
then it can express the functions Zi(t), (i = 2, . . . , r) through the z±(t) function. 
It may observe a relationship in eqn. ( |13| ) for z 2 (t) and Z\(t) functions respec- 
tively z 2 (t) = T L>™ 1_m2 ~ ai+Q2 zi(i). We assume an equalities of the integer 
orders mi = m 2 . The previous relationship becomes z 2 (t) = T D^ ai+a2 z\{t). 
Following the property ( T D^ a y){t) = ( T If ?/)(*) we have z t (t) = T / t Q1_Ql zi(tj, 
(i = 2, . . . , r). In the ordinary differential equation presented in the system (y_3|) 
we change all Zi(t), (i = 2, . . . , r) which depends on z\(t). Moreover, we introduce 
a new temporary function in the following form 

W (t) = Zl (t) + ^ T J^-^Zl{t) + ...+ — rl^-^Z^t). (19) 



We need to find an inverse form of eqn. (|T^) as 



Zl (t) = (1 + - T 7 t a 



0l 



')-^(t), 



(20) 



where (l + ^r/? 1 ' 
the operator (l + fa T /« 



St r / t a i «»-) 1 denotes the left inverse operator to 



.at jo 



'). We can apply the known Laplace 



1 H| M/i 

transform for eqn. (119) that to find the left inverse operator in expression 
The formula (EOT) is not suitable for practical purposes and computations. 



Proposition 3. Connecting expressions ( ]J 3( j and Qjlty^r ($0\), the solution of the 
second subclass of the fractional differential equation is presented as follows: 



f D^w{t) = -i(o r+ i D™ r+1 2 I . + i(i) + . . . 
+a„ £> t m »z„(t) + o„+i y{t)) + £/(t) 
Zl (t) = (l + ^ a ir- a2 +...+ 
z r+1 (t) = Dr i - mr+1 - ai+ar+1 z 1 (t) 



< 



£oir~ op )-v*) 



(21) 



z„(i) = oA Tni-TOB ~ ai+ttB «i(*) 

y (1) (t) = EJSi" 1 ^^ + oA mi ~ ai+1 «i(*) 

We do not need to proof the above proposition because it was done in the 
previous subclass of the multi-term class. 



3 Numerical treatment and examples of calculations 

In this section we propose an explicit numerical scheme. For numerical solution 
of the differential equation of integer order we apply one-step Euler's method || . 
There is no limits to extend above approach to the Runge-Kutta method. Ac- 
cording to results presented by Oldham and Spanier |Q we use numerical scheme 
for integral operator as 



1 t z i 



h a 



2r(l + a) 



Zo (i a - (i - id +zt+j2 a-, ((? + ir - u - 1) 



which is valid for arbitrary a > 0. We also use an algorithm given by Oldham 
and Spanier |Q for numerical differentiation. The algorithm depends on a range. 
For < a < 1 we have 



(1 - a)z 



i-l 
J=0 



Zi _ u+1) )((j + iy-<*-f-<*) 



(23) 



In literature |Q one may find the algorithm for highest values. 

For limited practical solution of cqn. ( pp| ) we found in literature Babenko's 
method Q. In general form of eqn. ( |20| ) we cannot consider the left inverse op- 
erator. For that, we found only a simplified formula limited to the two fractional 
orders cxi, 0/2, which binomial expansion we can show as 

00 

zu = T^(-l) j (-y • I {ai - a ^ j Wi. (24) 

3=0 

For some comparison we select the most popular fractional differential equa- 
tion in literature 

a -o D 2 t y{t) + b -o D\ A y(t) + c ■ y(t) = f(t) (25) 



called the Bagley-Torvik equation. Fig. |lja presents the qualitative comparison 
of the numerical solution generated by our method with the analytical solution. 




2 4 6 8 10 12 14 16 18 20 2 4 6 8 1 1 2 14 1 6 18 20 



Fig. 1. Solution of the Bagley-Torvik equation: 

a) comparison of our results with the analytical solution ^ and with Podlubny's 
work g, 

b) numerical solution of the non-linear form of Bagley-Torvik equation. 



initial conditions y(0) = 0, y'(0) = respectively. The analytical solution one 
may find in M. We also found a simple way of solution created in Podlubny's 
work g . Fig.|j]a certifies that our numerical results behave good in comparison to 
the analytical solution. The great advantage would be the apply of our approach 
to the non-linear form of fractional differential equation. Let us consider the 
non-linear form of the fractional differential equation 

Q D 2 t y{t) + 0.5 • oA 15 yW + 0-5 • y\t) = f(t) (26) 

where f(t) = q 1 < t < oo wn ^ initial conditions y(0) = 0, y'{0) = 0. This 
is the Bagley-Torvik equation where non- linear term y 3 (t) is introduced. Fig. |l|b 



shows a behaviour of the numerical solution of eqn. (pff). We can see that the 
step h has strong influence to the solution y(t). Non- linear fractional differential 
equations need some computational tests that to choose right value of the step. 

4 Conclusions 

In this paper we propose a new method for the numerical solution of arbitrary 
differential equations of fractional order. Following to Blank's results the main 
advantage of the method is a decomposition of fractional differential equation 
into a system composed with one ordinary differential equation of integer order 
and the left inverse equations of the Abel- integral operator. We distinguish two 
classes of such system: one-term fractional derivative and multi-term fractional 
derivatives. The comparison certifies that our method gives quite good results. 
Summarizing these results, we can say that the decomposition method in its 
general form gives a reasonable calculations, is the effective method and easy to 
use and is applied for the fractional differential equations in general form. 
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